Signatures of nonadiabatic O2 dissociation at Al(lll): 
First-principles fewest-switches study 



Christian Carbogno,^ Jorg Behler,^ Karsten Renter,'^ and Axel Grofi^ 

' Institut fiir Theoretische Chemie, Universitdt Ulm, D-89069 Ulm, Germany 
^Lehrstuhl fiir Theoretische Chemie, Ruhr- Universitdt Bochum, D-44'^80 Bochum, Germany 
^ Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, D-14195 Berlin, Germany 

(Dated: October 2, 2009) 

Recently, spin selection rules have been invoked to explain the discrepancy between measured 
and calculated adsorption probabilities of molecular oxygen reacting with Al(lll). In this work, 
we inspect the impact of nonadiabatic spin transitions on the dynamics of this system from first 
principles. For this purpose the motion on two distinct potential-energy surfaces associated to 
different spin configurations and possible transitions between them are inspected by means of the 
Fewest Switches algorithm. Within this framework we especially focus on the infiuence of such 
spin transitions on observables accessible to molecular beam experiments. On this basis we suggest 
experimental setups that can validate the occurrence of such transitions and discuss their feasibility. 
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I. INTRODUCTION 

Metal oxidation processes are of outstanding im- 
portance in various industrial fields which range from 
combustion catalysis to microelectronics. Accordingly, 
this research field has attracted considerable scientific 
interest — both from experimentalists and theorists — over 
the last decades. Although the most fundamental rules 
underlying such processes were understood^ quite early, 
there still exist some not yet elucidated phenomena, even 
for systems that appear utterly simple at first sight. One 
of the most prominent and intriguing examples is the 
oxidation of the (111) surface of aluminum which ex- 
hibits various astonishing features. The exposure of a 
clean aluminum surface to an oxygen atmosphere leads 
to the growth of an oxide film, which then inhibits fur- 
ther corrosion. In spite of the fact that the overall mech- 
anism for the formation of this AI2O3 stoichiometry has 
already been clarified by qualitative^ kinetic considera- 
tions, the elementary key steps underlying this process 
are still topic of a lively scientific debate. 

For instance, scanning tunneling microscopy (STM) 
images obtained after exposing the clean aluminum sur- 
face to a minute dosage of molecular oxygen show only 
single oxygen atoms in an average distance of more than 
80 A from each other'^ on the surface. Various hypothe- 
ses, e.g., the subsurface migration of single oxygen atoms^ 
or a cannonball like abstraction dynamics^i^ of the im- 
pinging oxygen molecules, have been proposed to explain 
this extraordinary finding. As a matter of fact, even the 
validity of the original experiment has been questioned^. 
The contradictory nature of the proposed mechanisms, 
which are all backed up by reasonable scientific findings, 
already shows how complicated the experimental clarifi- 
cation of this issue is. 

However, at first sight the questions should be settled 
easily by means of the current available simulation tech- 
niques^. Unfortunately, this is not the case for this par- 



ticular system: Various studies carried out on the ba- 
sis of state-of-the-art adiabatic density functional the- 
ory (DFT )^°d^ were not even able to qualitatively re- 
produce the low adsorption probability for thermal O2 
molecules at the clean aluminum surface found in molec- 
ular beam experiments^. If even the simulation of such 
an elementary key step of the reaction fails that dramat- 
ically, attempts to model the more complex features of 
the dissociation dynamics within the framework provided 
by DFT are foredoomed to fail, too. 

Recently, however, DFT calculations addressing poten- 
tial energy surfaces (PESs) of fixed spin states were able 
to reproduce the experimental data in a semiquantita- 
tive way. These studies suggest that spin selection rules, 
which hinder^- an adiabatic spin transition from the ini- 
tial O2 gas-phase triplet state to the singlet state of the 
adsorbed O atoms, are the rate determining factor in this 
process. For these investigations a constrained density 
functional theory approach^^ was employed to compute 
the PESs of O2 in different spin-configurations, and the 
lowered adsorption probability was then correctly repro- 
duced when restricting"'^^ the O2 molecules to move on 
the spin-triplet potential-energy surface only. 

Although these studies yield a more than qualita- 
tive correspondence with the experiment with respect 
to the initial sticking coefficient, they inherently incor- 
porate an approximation that is ultimately not satis- 
factory. By restricting the motion to the triplet PES 
alone, the impinging O2 molecules are not able to re- 
lax to the correct final singlet state of the adsorbed 
atoms while dissociating. Therefore we extend these pre- 
vious studies by considering the dissociation dynamics on 
multiple spin potential-energy surfaces, allowing transi- 
tions between them within the mixed quantum-classical 
Fewest Switches^ algorithm as proposed by TuUy, which 
has been successfully employed in the modeling of var- 
iou a— 1 ^ molecule-surface processes before. As detailed 
below, we can thereby reproduce the measured sticking 
coefficient and thus reconfirm the earlier more restricted 
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modeling. Additionally, our simulations allow to study 
the occurrence, the nature and the influence of the spin 
transitions in detail. On the basis of these calculations, 
we were already able to propose^^ measurements that can 
clearly prove the nonadiabaticity of the studied reaction. 
In this work we further study such potential experimental 
setups and especially focus on their feasibility. Therefore, 
we explicitly inspect how various uncertainties underly- 
ing the theoretical modeling, as for instance the frozen 
substrate approximation, affect the practical realization 
of such scattering experiments in the lab. 



The proposed measurements can verify the fundamen- 
tal assumption behind our approach that hindered spin 
transitions are indeed the rate determining factor in this 
process. This is of particular importance due to the fact 
that speculations22i2i have arisen which do not relate 
the low sticking coefficient of thermal oxygen molecules 
to the nonadiabatic dynamics. Rather, shortcomings 
of the state-of-the-art but still approximative exchange- 
correlation GGA functional employed in the DFT cal- 
culations are held responsible for the discrepancy found 
between the adiabatic simulations and the experiment. 
Such conjectures are hard to confute theoretically, at 
least as long as calculations with an improved exchange- 
correlation functional are numerical too costly to be per- 
formed accurately on a large scale. Certainly, this dic- 
tates also further research along these lines. However, 
the herein proposed experiments can clarify this severe 
doubt already now and thus establish a solid founding 
for further research on the complex phenomena associ- 
ated with the oxidation of aluminum surfaces. 



In this paper we will first review the ab initio meth- 
ods underlying the simulations in Sec. [ill Since the DFT 
techniques used for the calculation of the total energy of 
the system as well as the techniques employed to inter- 
polate this set of data to a smooth potential-energy sur- 
face have already been discussed in detailiii^ before, we 
will focus on the method used to calculate the dynamics 
on multiple potential-energy surfaces, i.e., TuUy's Fewest 
Switches^'' algorithm. Especially, we will discuss the na- 
ture of the simulated electronic states and the nonadi- 
abatic coupling between them. In Sec. IIIII we will then 
present the results of such simulations, whereby a spe- 
cial focus lies on the occurrence and on the role of the 
nonadiabatic transitions. In this context we will also 
present possible experimental approaches able to iden- 
tify the presence of the spin transitions, whereby we will 
also discuss the feasibility of the proposed experiments. 
We especially focus on the uncertainties underlying the 
employed theoretical approach and try to quantify them 
by additionally inspecting the influence of different cou- 
pling strengths and the role of the surface mobility. 



II. THEORETICAL METHODS 

A. Mixed Quantum-Classical Dynamics 

When inspecting the dynamics of a system that in- 
cludes nonadiabatic transitions, the electronic and nu- 
clear degrees of freedom cannot be tackled independently 
and subsequently, as done in the Born-Oppenheimer ap- 
proximation. Unfortunately, a full quantum dynamical 
description of both the electronic and nuclear degrees of 
freedom is computationally far too demanding for the 
present system. However, for the nuclei usually a classi- 
cal treatment is sufficient. Still, the fast electronic and 
slow nuclear degrees of freedom have to be inspected si- 
multaneously, whereby a mutual self-consistent feedback 
between them is crucial^^ to achieve a correct description 
of the system. In this work we employ a mixed quantum- 
classical method, namely the Fewest Switches Surface 
Hopping algorithm as proposed by TuUy^*', to achieve 
this goal. For a concise summary of this approach, the 
total Hamiltonian for a system including electronic (r) 
and nuclear (R) degrees of freedom 

H = Tn + T, + Vr + VrR + VR 

= TR + i/,i(r,R) + VR (1) 

will serve as a starting point. In Eq. ([T]), T denotes the 
kinetic energy operators, V the respective electrostatic 
potentials and all terms depending on the electronic co- 
ordinates are subsumed in Hei ■ The total wave function 
for this system can then be expressed as 

vI'(r,R,t) = ^c,(t)<i>,(r,R) , (2) 

i 

whereby the basis functions ^^(r, R) used in the expan- 
sion, also referred to as electronic states in this work, are 
not necessarily eigenfunctions of the electronic Hamilto- 
nian Hei- For each of these basis functions, we can derive 
a Hamilton- Jacobi equation for the nuclear dynamics, in 
which the potential 

y,,(R) = ($,(r,R)|Fe/(r,R) + VR|$,(r,R)) (3) 

resulting from the integration over the electronic degrees 
of freedom indicated by the bra-ket notation, is specific 
to the thereby employed basis function $i (r, R) . Conse- 
quently, an equation of motion for each particular elec- 
tronic state included in the expansion ([2]) can be derived. 
Whenever the nuclear masses of interest are significantly 
larger than the one of hydrogen, so that quantum effects 
in the dynamics of the slow degrees of freedom can be 
safely neglected^!, the motion can be described by clas- 
sical equations of motion 

M^R=-VrF,,(R). (4) 

Still, these equations do just describe the motion on 
one single PES but not the transitions between different 
PESs. 
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If we insert the full expansion ^ into the time- 
dependent Schrodinger equation associated with the 
Hamiltonian H, we gain a set of coupled differential equa- 
tionsi^iii 



ihck — 



5^ c,{t) {Vk,{R{t)) - zhKit) ■ d,,(R(i))) , (5) 



which describe how the expansion coefficients related to 
the individual electronic states evolve along a trajec- 
tory R(t) . Beside the non-diagonal potential matrix ele- 
ments Vfcj(R(t)), also the nonadiabatic coupling vector 



d,,(R) 



($,(r,R)|VR|$,(r,R)) 



(6) 



which designates the dependence of the chosen basis 
functions $i(r,R) on the nuclear positions, enters the 
equations. By integrating Eq. (O along a trajectory, 
the changes in the occupation of the different electronic 
states during the dynamics, i.e., the time evolution of 
the diagonal elements of the density matrix aij{t) = 
Ci{t)*Cj{t), can be inspected for a specific trajectory. 

This fact can be exploited to combine Eq. (j4|) and ([5|) 
to a surface hopping method: At each time step At, the 
system is in exactly one electronic state k that deter- 
mines which potential Vkk is used in Eq. (jH) to compute 
the classical trajectory. On the other hand, the com- 
plete set of equations ([5]) is integrated along the trajec- 
tory to determine the occupation in all electronic states 
of the expansion. If by any chance another electronic 
state from the one that actually determines the nuclear 
dynamics becomes dominant, the trajectory "switches" 
to this state: From then, this new state will determine 
the potential in the nuclear equation of motion, as long 
as no other switch occurs. To actually decide whether 
and when such a "surface hop" occurs, a Metropolis-^ 
like algorithm is combined with the switching probabil- 
ity proposed by TuUyi^ 



it) 



a-kk 



At 



(7) 



This so called Fewest Switches approach, which distin- 
guishes this method from comparable surface hopping 
approaches (see Ref. [23| and references therein), fulfills 
two essential conditions. For a large number of com- 
puted trajectories, the population in each state statis- 
tically matches the probabilities given by the density 
matrix an. Thereby, this distribution is achieved with 
as few surface hops as possible. Please note that this 
algorithm allows switches to occur at any point along 
the trajectory R(t), even if the potential energies dif- 
fer (e.g., AV = Vjj{R) - Vkk{R) ^ 0). For such sur- 
face hops with Ay 7^ 0, special care has to be taken 
with respect to energy conservation. If the kinetic en- 
ergy of the nuclei is smaller than the potential difference 
(Tr, < AV), the switches are rejected. Conversely, if the 
potential energy difference is larger than the kinetic en- 
ergy (Tr > AV), the velocities have to be rescaled to 
-AV. This constraint is not unique 



fulfill 



because the direction of the velocity adjustment can be 
chosen freely. Considerations from semi-classical theory 
suggest to rescale the velocities along the nonadiabatic 
coupling vector, since the forces associated with nonadi- 
abatic transitions typically^^'^^ point in its direction. 

Although the herein given description of the Fewest 
Switches method is adequate to understand the calcula- 
tions presented below, it is certainly not exhausting. The 
interested reader is referred to the original publications 
describing this metho d^^i^^'^^ . 



B. The inspected electronic states 

The first essential step in the simulation of nona- 
diabatic reactions with the surface hopping method is 
the identification of the electronic states of interest, 
i.e., the basis used for the expansion Since a distinct 
PES Vii{Il), which has usually to be calculated individ- 
ually under huge numerical effort, is associated to each 
one of these electronic states, it is desirable to achieve 
a correct description of the system with as few expan- 
sion terms as possible. Fortunately, the inclusion of the 
adiabatic ground state and one additional excited state 
is typically already sufficient to model nonadiabatic pro- 
cesses accurately^^. In our case an analysis of the adi- 
abatic calculation — a so called diahatization — serves as 
a starting point for the determination of the electronic 
states that are relevant in the inspected system's dy- 
namics, since such a diabatia^ representation allows for 
a more natural and intuitive interpretation of the occur- 
ring electronic transitions. 

As shown for an exemplary geometry in Fig.[Tl the PES 
associated with the initial state of the reaction, i.e., the 
triplet ground state of the isolated oxygen molecule, re- 
sembles the adiabatic ground state PES closely at large 
molecule-surface distances and is thus the first state in- 
cluded in the surface hopping approach. Not too sur- 
prisingly, this triplet PES does not coincide with the 
adiabatic PES in the limit of small oxygen-surface dis- 
tances: In this final state of the reaction, the dissoci- 
ated oxygen atoms adsorbed on the surface exhibit a 
singlet type ground state electronic structure, whereby 
an additional charge transfer-iSiil from the surface onto 
the molecule occurs. Although there might be multi- 
ple charged states of the oxygen molecule involvedS^i in 
the actual electron transfer process, we can qualitatively 
incorporate such effects in this schematic model by one 
effective PES. This so called ionic PES resembles the adi- 
abatic potential closely in the late part of the dynamics 
and approaches the lowest lying charge-transfer state in 
the vacuum limit. Its respective energy gap can be esti- 
mated from the clean aluminum surfaces work function 
and the isolated oxygen molecules electron affinity to be 
4 eV with respect to the triplet ground state. In between 
these two states, two additional singlet states of the neu- 
tral oxygen molecule, i.e., the ^ Ag and the state, are 
located. Although these two states do neither correspond 
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FIG. 1: Schematic plot of the potential-energy surfaces along the minimum energy path in the triplet state for dissociation 
in parallel alignment over the fee site: Beside the adiabatic (spin-polarized), the triplet (constrained spin-polarized) and the 
singlet (adiabatic non-spin-polarized) PESs, which have been determined by means of DFT, also the expected functional form 
of the potentials associated with the first two excited singlet states and the so-called "ionic" state are shown in Fig. (a). In 
Fig. (fe) the lower and the upper bound derived for the coupling of the triplet and the singlet PES are plotted along the minimum 
energy path of Fig. (a). 



to the initial nor to the final state of the reaction, they 
play a crucial role in the dynamics. In the proximity of 
the barriers present on the triplet PES they definitely are 
as relevant as the respective triplet and ionic PES, since 
they exhibit comparable curve crossings with the triplet 
state. 

It would certainly be desirable to include all of these 
excited states in the simulation of the nonadiabatic dy- 
namics. However, even the first principles calculation 
of one single PES alone is utterly costly in terms of 
computing power. But even more importantly, the ^Ag 
and the ^T,^ state exhibit a strongly correlated, multi- 
determinant wave functioni^ii^i^ that is not described 
correctly by single-determinant theories such as Hartree- 
Fock and DFT. As a matter of fact, the current exchange- 
correlation functionals, which do not depend on the sym- 
metry of the electronic structure, do not allow to discrim- 
inate^!^ between the ^ Ag and the state at all. This 
becomes evident already when inspecting the isolated 
oxygen molecule by an adiabatic non-spin-polarized DFT 
calculation, which must yield a singlet state by construc- 
tion. In this approach, however, the respective wave 
functions of the two singlet states mi^s^ in equal parts, 
so that the resulting triplet-singlet gap lies between the 
correct values that are associated with the ^Ag and the 
-'^E+ state, respectively. When additionally including the 
aluminum slab in such an adiabatic non-spin-polarized 
DFT calculation, the actual triplet-singlet gap is lowered 
due to a spurious and unphysical charge transfer from the 
surface that occurs even at large molecule-surface separa- 
tions: The current GGA exchange-correlation functionals 
erroneously favor such a charge separation due to their 
incorrect convexity with respect to fractional charges"^^. 
As a consequence the thereby computed PES mimics the 
first excited ^Ag singlet state at large distances from the 



surface, as shown in Fig. [T] In the proximity of the sur- 
face, this PES coincides with the one associated with the 
adiabatic ground state, since the charge transfer between 
oxygen and aluminum is not restrained in the underly- 
ing adiabatic non-spin-polarized DFT calculations. Thus 
such an approach enables us to elfectively describe both 
the charge transfer effects close to the surface as well 
as the influence of the neutral singlet states within one 
single PES. However, such an approximative but still ad- 
equate description of the manifold of excited states also 
has some drawbacks. For the calculation of the coupling 
elements we have to disentangle the contributions from 
the various electronic states once again, as discussed in 
the next section. 

For the sake of readability the adiabatic spin-polarized 
PES will henceforth be referred to as "adiabatic", the 
constrained spin-polarized triplet PES as "triplet" and 
the adiabatic non-spin-polarized PES as "singlet". All 
three potential energy surfaces, which explicitly depend 
on the six degrees of freedom of the oxygen molecule, 
have been constructed in the frozen substrate approxi- 
mation by interpolating large sets of ab initio data. A 
detailed discussion of the underlying constrained and un- 
constrained DFT simulationsi^ii^, of the neural network 
interpolation^^ and of the dynamics on a single PESi^ii^ 
can be found in the cited publications. 

C. The Coupling 

For the derivation of the coupling elements we will 
again rely on the diabatic representation (see Sec. Ill B[) 
of the triplet and the singlet state (|^'t) , l^'s)), for which 
the nonadiabatic coupling vector 

dt,(R) = (^-tlVRl*,) =0 (8) 
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vanishes exactly at each point in space due to the trivially 
fulfilled orthogonality of the respective spinors. Conse- 
quently, the evolution of the density matrix in Eq. ([5]) — 
and hence the electronic transitions — is completely de- 
termined by the non-diagonal elements of the electronic 
Hamiltonian 



Vts - (^-tlHel*. 



(9) 



To disentangle the various interactions contributing to 
this coupling, we decompose the singlet wave function, 
which is actually a mixture of three distinct states (see 
Sec. Ill Bj) . into its individual components: 



A' 



V2 



B ■ 



(10) 



i.e., the first excited ^Ag singlet state the second 
excited singlet state VPs and the "ionic" state ^'io„, 
which adiabatically incorporates all charge transfer 
effects and is identical to the singlet ground state close 
to the surface. Thereby the complex parameters A and 
B are determined by the adiabatic mixing of the neutral 
and the ionic state described before. 

In a first step, we neglect the contributions arising from 
the coupling to the ionic state and focus on the isolated 
molecule in the vacuum limit, where no charge transfer 
is supposed to happen. Under these assumptions the 
coupling can be simplified to: 



1 

71 



((*t|We|*A) + (*t|He|*s)) 



(11) 



Since the spinors present on the left and on the right 
hand side of the Hamiltonians in Eq. (|lip are orthogonal 
to each other, the only interaction that does not vanish 
in the evaluation of this matrix element is the spin-orbit 
couphng (SOC) 

Vts = ^ ((*t| Vsoc I^-a) + (*t| 14oc |*s)) . (12) 



The first of the two coupling terms in Eq. (|T2l) . i.e., the 
SO coupling of the triplet ground state to the first ex- 
cited ^ Ag singlet state, vanishes in first order approxima- 
tion, since selection rules prohibit a direct mixing of the 
two states^i^^. Solely, an indirect couplin g^^i^^ through 
the highly excited '^Ilg state, which is at least two or- 
ders of magnitude smaller than the direct SO coupling 
of the triplet ground state to the "'^E+ state discussed 
below, is possible. Thus the overwhelming portion of 
the spin-orbit coupling is provided by the second term, 
i.e., the SO coupling of the triplet ground state to the 
second excited singlet state. For the isolated oxy- 
gen molecule, this non-diagonal matrix element has been 
calculated as a function of the oxygen-oxygen distance 
with the SOC module^S provided by the quantum chem- 
ical simulation package Molprc^'^ on the basis of wave 
functions determined^ by MuUireference Configuration 



Interaction'^^ before. The thereby determined coupling 
Vts — as well as the derived lifetime of 11.16 s for the 
state — is in excellent agreement with experimental and 
theoretical data found in literature^i^^'^. With respect 
to the oxygen-oxygen separation, which is the only de- 
gree of freedom the coupling actually depends on in these 
calculations, it is a linear, slowly decreasing (almost con- 
stant) function, even up to huge separations (see Fig. [1]). 
If actually a factual dissociation of the isolated molecule 
does not affect this coupling, also the interaction with 
the aluminum, which exhibits only a small SOC— due to 
its low mass, will not alter it significantly. Consequently, 



(13) 



the absolute value of which is plotted in Fig. [1] can be 
regarded as a lower bound for the true coupling present 
in the system. 

On the basis of this minimal coupling determined 
above and the PESs described in Sec. Ill Bl we can con- 
struct the diabatic potential matrix 



Et K„.„ 



V* 



Es 



(14) 



and diagonalize it to obtain the respective adiabatic rep- 
resentation of the potentials 



E,. 



Et±^JiEs- EtY +4.\V^.^\ 



(15) 



for this minimal coupling. Due to its minute value, how- 
ever, only a small deviation from the diabatic potentials 
occurs at all, as already evident from the minute gap of 
approximately 2t^„i„ w 30 meV at the avoided crossing 
point. This fact also highlights the main shortcoming of 
this minimal coupling element: The minute splitting at 
the avoided crossing is not able to recover the adiabatic 
ground state potential that has been determined by the 
unconstrained DFT calculations before. However, this 
fact also allows to turn the tide: By forcing the left hand 
side of Eq. to correspond to the adiabatic potential 
determined by DFT before, we can invert the diagonal- 
ization to obtain a complementary value for the coupling 
element 



K,.x| = \EsEt~EsEl 



EtE°r. + {E^^jf. (16) 



In contrast to the one-dimensional minimal coupling, this 
coupling term |K,iax| is a true six-dimensional function 
of the oxygen coordinates, since the inverse diagonaliza- 
tion (|16p can be evaluated at each point in space for which 
respective triplet, singlet and adiabatic potential values 
are available. Accordingly, this coupling also depends on 
the respective rotational and center of mass degrees of 
freedom so that it does not exhibit such a constant be- 
havior along the minimum energy path as the SOC (see 
Fig. [T]). Due to the spurious charge transfer occurring 
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in the adiabatic calculations, however, this maximal cou- 
pling does not vanish in the vacuum limit as one would 
expect. Additionally, the inverse diagonalization proce- 
dure itself is generally expected to overestimate the cou- 
pling^*^, so that Ki,ax can be regarded as an upper bound 
for the true coupling present in the system. In fact we 
believe that the maximal coupling is closer to the true 
coupling strength than the minimal coupling. 

In the case of such a maximal coupling, a diagonaliza- 
tion of the respective diabatic Hamiltonian matrix yields 
a ground state potential that resembles the adiabatic 
PES exactly, as required by construction. Consequently 
all barriers previously present in the diabatic picture van- 
ish as well, so that surface hopping simulations of ther- 
mal molecules, for which the diabatic calculations yield a 
sticking probability of 2% (see Sec. IIII Ap . lead to a con- 
stant sticking coefficient of 100% in the adiabatic picture. 
This alarming discrepancy between the two representa- 
tions can be traced back to our approximative description 
of the manifold of excited states by a single effective PES. 
In the limit of a large coupling our model describes the 
two-level system consisting of the triplet and the ionic 
state, but does not account for the presence of the other 
two neutral singlet states. 

In the diabatic picture this is indeed a justified ap- 
proximation, since due to the minute coupling the two 
singlet states just provide for a negligible amount of addi- 
tional nonadiabatic transitions with respect to the ionic 
state. This is no longer true in the adiabatic picture: 
In this representation, the potential-energy surfaces of 
states that are strongly coupled are also strongly dis- 
torted at avoided curve crossings. The smaller the cou- 
pling of the respective states is in the diabatic picture, 
however, the more probable transitions between the adi- 
abatic PES become. By not including the singlet states 
and their minute coupling in the adiabatic representa- 
tion of our two-level model, we are de facto ignoring the 
adiabatic PES that matters most, i.e., the one to which 
a transition will almost certainly occur at the avoided 
crossing. To avoid these problems without the cumber- 
some construction of another PES (see Sec. Ill B( ). we per- 
formed all the calculations presented here in the diabatic 
representation. The price to pay is that the effective bar- 
riers in the simulations are slightly overestimated, since 
nuclear tunneling effects and resonances, which would 
be inherently incorporated^ in an adiabatic picture, are 
neglected in such a diabatic representation. Fortunately, 
these effects are typically negligible in high-dimcnsional, 
heavy systems^i^ such as the one addressed in this work. 



D. Calculational Details 

In all presented MD and SH simulations the initial 
conditions are chosen to model a molecule in the vac- 
uum limit, in which it does not interact with the alu- 
minum surface. Therefore the initial center-of-mass dis- 
tance between the oxygen molecule and the surface is 



always 5 A. Also, a certain amount of energy E is as- 
signed to the translational motion perpendicular to the 
surface. The initial oxygen-oxygen distance and vibra- 
tion is chosen to model the zero point energy of the free 
molecule (« 0.1 eV) in the so called quasiclassical ap- 
proximation^i. In the beginning no kinetic energy is 
assigned to the rotational and to the lateral center-of- 
mass degrees of freedom, but their initial values are also 
chosen randomly to uniformly examine all possible ini- 
tial conditions. To sample the whole phase space and 
to converge the calculations with respect to the elec- 
tronic transitions at least 2, 000 trajectories have been 
computed for each of the data points in the following fig- 
ures. When a higher accuracy was needed, as in Sec. lIIIBl 
for instance, up to 2, 000, 000 trajectories have been com- 
puted. In each case the calculation of the individual tra- 
jectories is stopped when one of the following conditions 
is reached: Molecules are regarded as dissociated, if the 
oxygen-oxygen distance exceeds 2.0 A 160% of the 
equilibrium bond length of the free molecule), and as 
reflected, if the center-of-mass distance to the surface ex- 
ceeds 5 A. The maximal time interval for one trajectory 
is chosen so that virtually all molecules eventually fulfill 
one of these two conditions. 

Naturally, the convergence of the results with respect 
to the presented calculational details and with respect 
to the numerical integration routines has been examined 
extensively^. 



III. RESULTS AND DISCUSSION 

A. The Sticking Coefficient 

The sticking coefficient So{E) at normal incidence, 
i.e., the probability of dissociative adsorption as a func- 
tion of the incident translational energy E for a molecular 
beam oriented perpendicular to the surface, is the basic 
observable that we are going to discuss while investigat- 
ing the nonadiabatic dynamics of this system. As shown 
in Fig. [21 it exhibits an S-like shape in the experiments — 
a typical evidence for an activated process. MD cal- 
culations on the adiabatic PES alone, however, yield a 
sticking coefficient of constantly 100% even at thermal 
incident energiesi^ due the almost complete absence of 
barrierai^i^ on the respective PES. Conversely, MD sim- 
ulations on the triplet PES reproduce the measurements 
more than satisfactorily, as already discussed in the in- 
troduction. 

With respect to the MD simulations performed on the 
triplet PES alon o^'^i^^ , the inclusion of nonadiabatic tran- 
sitions within the surface hopping method changes the 
simulated sticking coefficient only quantitatively, but not 
qualitatively. At large incidence energies there is ob- 
viously no substantial difference at all: The impinging 
oxygen is able to overcome the barriers anyway, regard- 
less of any occurring electronic transitions. In contrast to 
this, molecules at medium velocities are not always able 
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FIG. 2: (a) Sticking coefficient So{E) as computed by moiecuiar dynamics (MD) on tlie tripiet potentiai-energy surface and 
by Tuily's surface flopping metliod (SH) for both a maximai and a minimal coupiing of tlie tripiet and tlie singlet PES. Tlie 
respective experimental data^^ is shown as weli. (6) Area-normalized histograms of the barriers found on the individual PESs 
by the two-dimensional Nudged Elastic Band technique (see text). The potentials along the minimum energy path over the 
bridge site are shown in the inset. 



to overcome the barriers on the triplet PES. Whereas 
such molecules are repelled in the triplet MD calcula- 
tions, they can dissociate by switching to the singlet PES 
in the SH simulations, as soon as they reach the respec- 
tive crossing point (see Fig. [T^) . The probabihty for such 
an electronic transition to occur depends on the coupling 
of the two electronic states so that the relative increase in 
dissociation found in the SH simulations varies with the 
coupling strength. Last but not least, the triplet MD and 
the SH simulations yield the same sticking coefficient at 
the smallest incident energies below 150 meV: Molecules 
that are not even capable of reaching the crossing point 
are also not able to benefit from the possibility of an 
electronic transition. This might seem surprising at first, 
since one may expect at least some triplet-singlet crossing 
points to occur also in this low energy regime. However, 
this is not the case: As exemplarily shown for the bridge 
site in the inset of Fig.l^b, the crossing point and the top 
of the barrier almost coincide even in geometries that ex- 
hibit barriers well below 150 meV. As a consequence, the 
sticking coefficient is not affected by the switches in these 
setups, in spite of the fact that there is a high probability 
for such transitions at these low velocities. 

These effects can be studied more systematically on 
the basis of the actual barrier heights that the molecules 
experience on the individual PESs in the sudden approx- 
imation. To determine this statistical distribution we 
inspected the minimum energy path on a large number 
of elbow plots by a two-dimensional, self-implemented^^ 
Climbing- Image Nudged Elastic Band^^'"^^ method. Both 
the histograms of the barriers found on the adiabatic PES 
and on the triplet PES are shown in Fig. [2]d. The effec- 
tive barriers experienced by the impinging molecules in 
the surface hopping simulations, however, depend on the 
probability for an electronic transition and thus both on 



the coupling and on the velocity of the molecules. At fast 
velocities and for a minute coupling, the molecules stay 
on the triplet PES and thus behave as in the triplet MD 
simulations. In this limit, the barrier distribution of the 
triplet state describes the dynamics best. Conversely, the 
impinging molecules switch to the singlet state almost in- 
stantaneously while bypassing the crossing seam at low 
velocities and for huge couplings. In this limit, the ener- 
getic location of the crossing points can be regarded as 
an effective barrier distribution. 

As shown in Fig. [2b, the average height of these ef- 
fective barriers has also been determinedSS. by means of 
the Climbing Image Nudged Elastic Band method. The 
resulting barrier distributions substantiate the previous 
findings of the simulations: The adiabatic PES is unac- 
tivated for the large majority of geometries and even in 
the case of a hindered dissociation the arising barriers 
are typically smaller than 0.5 eV. In contrast to this, the 
triplet PES exhibits exclusively activated entrance chan- 
nels with a broad distribution of barriers that ranges 
from approximately 40 meV up to more than 1.5 eV. 
These barriers are drastically reduced by the occurrence 
of electronic transitions at the crossing points, the en- 
ergetic height of which is 0.5 eV at most. However, the 
lowest occurring barriers are not modified due to the rea- 
sons discussed before. 

Please note that the great majority of these barri- 
ers lies in the entrance channel, i.e., before the oxygen 
molecule starts to dissociate and the oxygen atoms move 
apart. Accordingly, no distinct influence of the vibra- 
tional state on the dissociation probability is expected. 
Nevertheless, the slight vibrational enhancement found 
in the experiments^ at small incident energies can be re- 
produced^ by the simulations and can be traced back to 
an effective lowering of the barriers due to an adiabatic 
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E(eV) 

FIG. 3: Reflection coeflicient at high incident translational 
energies as computed by MD on the triplet PES and by Tully's 
surface hopping method (SH) for a minimal and a maximal 
coupling of the triplet and the singlet PES. 



energy transfer— occurring in tlic nuclear dynamics from 
the vibrational degree of freedom onto the translational 
one. 



B. High Energy Scattering 

Although the hitherto presented simulations reproduce 
the experimental data successfully, this does not ex- 
plicitly prove the hypothesis^^ that nonadiabatic spin- 
transitions indeed play a decisive role in the dynamics 
of the triplet oxygen molecules at the aluminum sur- 
face. Therefore we searched for a possible experimental 
setup that is able to unambiguously demonstrate the oc- 
currence of such transitions. The fact that the triplet 
PES exhibits barriers as high as 1.5 eV might indeed 
offer a route towards experimental verification: For inci- 
dent energies that are larger than the triplet-singlet gap 
(« 1 eV), energy conservation does no longer prohibit 
a conversion of triplet into singlet molecules during the 
backscattering. If such singlet molecules were detected 
in the reflection channel, the nonadiabaticity of the scat- 
tering process could be experimentally proved. 

To investigate the feasibility of such a measurement, 
we focus on the reflection coefficient Ro{E) = 1 — So{E), 
i.e., the probability for a molecule with perpendicular in- 
cidence to be repelled back into the vacuum as a function 
of the initial translational energy E. As shown in Fig. [31 
we indeed find a small amount of molecules (< 5%) that 
are reflected in the triplet MD simulations due to the bar- 
riers on the triplet PES. The larger the incident energies, 
the less probable such backscattering events obviously 
become. For incident energies above 1.5 eV, for which 
the triplet PES does not exhibit any barriers anymore, 
the reflection coefficient does however not drop to zero 
but levels off at a value of 1.4 %. In this energy regime, 
the translational movement is so fast that the oxygen 



molecules approaching the surface in an unfavorable ge- 
ometry are not able to dissociate in the small time inter- 
val before being backscattered, since they do not follow 
the minimum energy path anymore. 

Not too surprisingly, the SH simulations with a min- 
imal coupling closely resemble the triplet MD calcula- 
tions. For such huge velocities and a minute coupling, 
oxygen behaves diabatically, i.e., it just stays on the 
triplet PES in almost all trajectories. In contrast to 
this, in the simulations with maximal coupling a signifi- 
cant amount of molecules still relaxes to the singlet PES, 
which exhibits no barriers, as soon as the crossing point 
is passed. Accordingly, the associated SH simulations 
yield only a minute reflection coefficient that slightly in- 
creases with larger incident energy. Even for such a 
huge coupling, the probability for electronic transitions 
decreases with larger velocities, so that for large incident 
energies more and more molecules stay on the more repul- 
sive triplet PES and are backscattered. Still, this effect is 
almost negligible, even for the fastest molecules. Regard- 
less of the coupling, the already small reflection proba- 
bilities are even further diminished as soon as one gives 
up the frozen substrate approximation (see Sec. IIII D[) . 
Even worse, less than 3% of these reflected molecules are 
found in the singlet state upon reflection, so that only 
0.01% of all incident molecules result in a backscattered 
singlet molecule at most, even in in the frozen substrate 
approximation. This fact makes an experimental detec- 
tion of the nonadiabatic transition in the presented way 
factually impossible. 



C. Singlet Scattering 

The sticking coefficient 5*0 {E) of triplet molecules dis- 
cussed above does not provide an unambiguous experi- 
mental evidence for the occurrence of nonadiabatic spin- 
transitions. This is principally caused by the fact that 
this typically inspected observable is not particularly sen- 
sitive to the occurring electronic transitions. For the 
same reason also a pure triplet MD approach yields a 
more than satisfactory agreement with the experimental 
data (see Sec. IIII At . To actually nail down the influence 
of the spin-transitions on the reaction dynamics exper- 
imentally, a complementary process that is significantly 
more sensitive to the transitions has to be found. In- 
deed, we could identifyiS, the dissociative adsorption of 
singlet oxygen molecules to provide for experimental ob- 
servables that strongly depend on the occurrence of such 
spin-transitions. A comparable strategy has already been 
used in a recent combined experimental and theoretical 
study^^ to demonstrate the importance of spin selection 
rules for O2 interacting with small Al anion clusters (~10 
to 20 atoms). In this work, an odd/even pattern in the 
reactivity of triplet oxygen has been found as a function 
of the number n of atoms in the Al clusters, whereas no 
such pattern has been found for the reaction of singlet 
molecules with these clusters. These studies motivated 
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FIG. 4: (a) Reflection coeflicient Rq as function of the incident 
kinetic energy of singlet oxygen molecules, as calculated by 
the SH method with a minimal coupling. Additionally, the 
respective data obtained by MD simulations on the singlet 
PES alone is shown. (6) Relative triplet yield Rt/Ro for 
both calculations. 



further studies addressing the dynamics of singlet oxy- 
gen, for example in the combination with other metal 
surfaces such as Ag(lOO)^. 

In our simulations we studied the dynamics of singlet 
oxygen molecules impinging in normal incidence on the 
aluminum surface at various kinetic energies. In the fol- 
lowing, we will concentrate on two complementary ob- 
servables that will allow us to understand the nonadia- 
batic processes active in the reaction. Firstly, we will 
inspect the total amount of backscattered molecules as 
function of the incident translational energy, i.e., the 
reflection coefficient Rq{E). In a second step, we will 
also inspect the relative amount of molecules that are 
backscattered in the triplet state, i.e, the relative triplet 
yield Rt{E)/Rq{E), which turns out as the actual fin- 
gerprint of the nonadiabatic transition. 

As shown in Fig. IH the probability of backscattering 
rises with increasing kinetic energy in the case of a mini- 
mal coupling: Although the singlet PES does not exhibit 
any barriers at all, some molecules at high velocities are 
reflected on this PES before the oxygen molecules are 
able to dissociate. The higher the incident energy, the 
more probable such a process obviously becomes. There- 
fore the respective SH simulations resemble MD simu- 
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lations on the singlet PES closely, apart from a con- 
stant offset. This difference is caused by the fact that 
in the backscattering process a reallocation of the trans- 
lational energy onto the other degrees of freedom occurs. 
Such vibrational or rotational excited molecules might 
not have the translational energy to escape the surface 
on the singlet PES anymore. While such molecules be- 
come trapped close to the surface and then finally disso- 
ciate in the singlet MD calculations, such molecules can 
escape in the SH simulations by switching to the triplet 
PES thus gaining the energy of the singlet-triplet split- 
ting which helps them to leave the surface. This mech- 
anism leads to an exponentially increasing reflection co- 
efficient at small incident energies, for which almost all 
backscattered molecules leave the surface in the triplet 
state. 

Intuitively one might expect the same mechanism to 
be active also in the case of the maximal coupling (see 
Fig. [5]). This is indeed the case at large incident energies, 
for which a small amount of molecules is backscattered 
on the singlet PES and then partially escapes the surface 
in the triplet state. The increased coupling favors the 
transition to the triplet state and thus increases both the 
total amount of reflected molecules as well as the relative 
fraction of triplet oxygen. A totally different mechanism, 
which is not active in the simulations with the minimal 
coupling, causes an exponential increase of the reflection 
coefficient for low and medium incident energies in the 
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simulations with the maximal coupling. When approach- 
ing the surface, the strong coupling destabilizes the sin- 
glet state, so that a decay to the triplet ground state be- 
comes extremely probable even while the crossing seam 
is approached. Thus the great majority of the molecules 
hit the surface on the triplet PES, on which there are 
barriers as high as 1.5 eV present (see Fig. [Sb), so that a 
backscattering event is orders of magnitude more likely 
than on the unactivated singlet PES. Consequently, one 
would expect the singlet reflection probability to be sim- 
ilar to the one of triplet molecules that hit the surface 
with a comparable kinetic energy E + AEts- 

A comparison with such a shifted triplet MD reflection 
coefficient surprisingly reveals that the respective reflec- 
tion probability of the singlet molecules is even higher, 
as shown in Fig. [5l During the nonadiabatic relaxation 
to the triplet state, which does not exclusively occur at 
the crossing point but also in its proximity where the 
two PESs do not coincide, the surplus of potential en- 
ergy is reallocated onto all degrees of freedom. Since this 
allotment is performed along the nonadiabatic coupling 
vector (see Sec.|ll|, not all the potential energy stored in 
the triplet-singlet gap E^^ is transferred onto the transla- 
tional degrees of freedom, which are principally respon- 
sible for surmounting the barriers. Rather, also the vi- 
brational, rotational and lateral center-of-mass degrees 
of freedom, which do not promote the dissociation as 
strongly or might even hinder it, gain a certain portion 
of kinetic energy, which is then missing in the perpen- 
dicular translational motion. Consequently, these sim- 
ulations are not directly comparable to the triplet MD, 
in which the complete amount E + AEt^ is assigned to 
the translational degree of freedom perpendicular to the 
surface. 

The fact that there are two different mechanisms ac- 
tive at low incident energies depending on whether the 
minimal or the maximal coupling is employed becomes 
evident when inspecting the velocities that the triplet 



molecules exhibit in the reflection channel. Fig. [6] shows 
the average kinetic energy gain that an incident sin- 
glet molecule experiences when being backscattered as 
a triplet. Since this data has been derived from the ve- 
locities that the reflected triplet molecules exhibit back 
in vacuum, the plots include both the energy rearrange- 
ment due to the nonadiabatic transition, but also all adi- 
abatic energy redistribution processes that occur during 
the backscattering dynamics. The two different mecha- 
nisms active at low incident energies result in a qualita- 
tively different reallocation of the energy. At an incident 
energy of 0.2 eV for instance, the strongest energy gain 
occurs in the translational degrees of freedom for the min- 
imal coupling, whereas it occurs in the rotational channel 
for the maximal coupling. 

The fact that the triplet-singlet conversion results in a 
kinetic energy gain also facilitates the realization of the 
proposed experiment. At low incident energies, there is 
no need to detect the electronic character of the backscat- 
tered molecules by, e.g., phosphorescence, since energy 
conservation allows to disentangle the contributions of 
the two molecular species by measuring the kinetic en- 
ergy distribution of the reflected moleculesS^, even in the 
case of a mixed triplet /singlet incident molecular beam. 
If any of the degrees of freedom exhibits an energy that 
is not compatible with the translational energy that the 
molecular beam initially exhibits, a transition from the 
singlet state must have happened. It is thus sufficient 
to perform a TOF measurement to unambiguously de- 
tect a significant amount of triplet oxygen through the 
respective kinetic energy gain. 



D. The frozen substrate approximation 

To further substantiate our findings, we also inspected 
the role of the surface mobility, which has not been ac- 
counted for in the hitherto presented calculations. As 
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FIG. 7: Reflection coefiicient Rq as function of the incident kinetic energy of singlet oxygen molecules, as calculated by the 
SH method with a (a) minimal coupling; (6) maximal coupling. The results both within the frozen substrate approximation 
(dashed lines) and for a three-dimensional surface oscillator (squares and triangles) are shown. In the respective insets, the 
relative triplet yield Rt/Ro for both calculations is plotted as well. 



mentioned before, the employed potential-energy sur- 
facesi^ solely depend upon the degrees of freedom of the 
oxygen but not on the ones of the aluminum atoms. In 
this so called frozen substrate approximation no energy 
dissipation into the bulk is possible, so that our descrip- 
tion of the dynamics close to the surface, where the ad- 
sorption energy of the oxygen molecule is freed, is cer- 
tainly defective. For the actual triplet dissociation prob- 
ability at low incident energies, which is basically deter- 
mined by the barriers present in the entrance channel (see 
Sec. lIII A| . the dynamics of the bulk atoms can generally 
be neglected. Still, the frozen substrate approximation 
might affect the reflection coefficients for the singlet oxy- 
gen, in the dynamics of which backscattering processes 
close to the surface play an important role. 

To overcome the frozen substrate approximation we 
performed simulations with an additional surface oscil- 
latori^j^^^, which can account for the possible energy 
transfer between molecule and bulk in hindsight, i.e., on 
the basis of potential-energy surfaces derived within the 
frozen substrate approximation. For this, the aluminum 
slab is permitted to oscillate as a whole in a harmonic 
potential, so that the actual molecule-surface distance 
does then not solely depend on the cartesian center-of- 
mass coordinate z of the oxygen anymore but also on the 
amplitude of the surface oscillator z. Thus both the di- 
agonal and the non-diagonal potential matrix elements 
that enter the surface hopping simulations in Eq. ^ and 
([5]) have to be formally adjusted to reflect this additional 
dependence: 



matrix 



Surf. Osci. 



■)• (17) 



Furthermore, the auxiliary harmonic potential of the 
slab, which does not depend on the electronic state, is 
added to each of the diagonal elements of the potential 



(18) 



whereas no additional coupling term is added to the non- 
diagonal elements. Analogously, this model can be ex- 
tended to also allow for an oscillation of the slab with 
respect to the lateral surface coordinates, as done in 
all the simulations presented below. Although such a 
surface oscillator does not yield any additional informa- 
tion about the actual motion of the individual aluminum 
atoms, it effectively incorporates the energy transfer ef- 
fects in between the slab and the molecule into the sim- 
ulations. For the frequency uj of the harmonic potential 
in Eq. (fT8|) a value of 70% of the experimental Debye 
frequency (« 34 meV) of the bulk^ has been used, since 
previous studies showed^^i^iS that this fraction is best 
suited to model the energy dissipation. For the exact 
same reason three times the mass of a single aluminum 
atom is assigned to the surface oscillator to account for 
finite size effects^i^. In the herein presented calcula- 
tions, a kinetic energy that corresponds to a temperature 
of 330 K has been assigned to the surface vibration in the 
initial condition. However, no significant dependence on 
this parameter has been found in the temperature regime 
at which the experiments are typically performed^i^ i^^'^" . 

Not too surprisingly, the sticking coefficient for inci- 
dent triplet molecules is hardly affected by the surface 
oscillator at all due to the fact that this observable is 
principally determined by the dynamics in the entrance 
channel. As a matter of fact, absolute changes of ±3% 
due to the inclusion of the surface oscillator have been 
found in the simulations at most^ and are thus not fur- 
ther detailed in this work. In the high energy scatter- 
ing regime discussed in Sec. IIIIBl however, the already 
minute reflection coefficients are even further diminished 
due to the surface mobility. 
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In contrast to this, the surface oscillator affects the re- 
flection coefficient for singlet molecules in a much more 
peculiar fashion that sheds light on the underlying tran- 
sition mechanisms: In the case of a maximal coupling the 
effects are again minute, as shown in Fig. [71 At large in- 
cident velocities, the surface mobility favors dissociation 
as intuitively expected, since the motion of the aluminum 
atoms offers a dissipation channel for the incident energy. 
At low velocities, the inclusion of the surface oscillator 
leads to a slight increase in the reflectivity. Already when 
climbing the barriers the molecules lose a small portion 
of their kinetic energy to the surface oscillator due to 
Newton's second law. In certain cases, this minute en- 
ergy transfer is already sufficient to hinder the surmount- 
ing of the respective barriers. However, the net effect is 
again minute, since these processes occur far away from 
the surface. 

We find a more pronounced influence of the surface 
mobility in the surface hopping calculations with a min- 
imal coupling, as also shown in Fig. [71 At large incident 
energies the reflection probability drops to zero. The in- 
cident molecules become trapped close to the surface due 
to the energy loss to the surface oscillator when hitting 
the aluminum, and thus eventually dissociate. The ex- 
act same mechanism is active at low incident energies, 
as also discussed before. However, the energy dissipation 
through the bulk slows down the oxygen even more, so 
that a transition to the triplet state becomes more and 
more probable in spite of the minute coupling value. Ac- 
cordingly, we find an increase in both the total and the 
relative yield when the mobility of the surface atoms is 
accounted for. 



IV. CONCLUSION 

In this work we extended previous MD studies^ of the 
dissociative adsorption dynamics of oxygen on the alu- 



minum surface by including multiple PESs associated to 
different spin-configurations in the dynamics by means 
of TuUy's Fewest Switches algorithm. Apart from re- 
producing the experimental sticking coefficient at normal 
incidence, this approach additionally allows to study the 
characteristics of the underlying nonadiabatic transitions 
in detail. This knowledge allows to propose experiments 
that are able to unambiguously identify the occurrence 
of such spin transitions. Whereas the detection of the 
triplet-to-singlet conversion in high energy scattering ex- 
periments is highly unlikely, a non-zero probability for 
the singlet-to-triplet conversion in scattering experiments 
with a beam of singlet oxygen molecules has been found, 
even in the limit of an unrealistically small coupling. If 
we additionally include the motion of the substrate atoms 
in our model, the determined yields are either not altered 
significantly or even enhanced. Depending on the cou- 
pling strength, two different mechanisms, which result 
in different kinetic energy distributions in the reflection 
channel, have been found in the singlet-to-triplet conver- 
sion dynamics. An experimental determination of these 
distributions does thus allow (a) to prove the nonadia- 
baticity of this reaction and (b) to draw conclusions on 
the character of the ongoing transitions. 
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